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Abstract 

The excitation energies and ionization potentials of the atoms in 
the first transition series are notoriously difficult to compute accu- 
rately. Errors in calculated excitation energies can range from 1-4 
eV at the Hartree-Fock level, and errors as high as 1.5eV are en- 
countered for ionization energies. In the current work we present and 
discuss the results of a systematic study of the first transition series us- 
ing a spin-restricted Kohn-Sham density-functional method with the 
gradient-corrected functionals of Becke and Lee, Yang and Parr. Ion- 
ization energies are observed to be in good agreement with experiment, 
with a mean absolute error of approximately O.lSeV; these results are 
comparable to the most accurate calculations to date, the Quadratic 
Configuration Interaction (QCISD(T)) calculations of Raghavachari 
and Trucks. Excitation energies are calculated with a mean error of 
approximately 0.5eV, compared with ~ leV for the local density ap- 
proximation and O.leV for QCISD(T). These gradient-corrected func- 
tionals appear to offer an attractive compromise between accuracy and 
computational effort. 







1 Introduction 



The accurate computation of the excitation energies and ionization potentials 
of the first transition metal series has proven to be a difficult problem for 
electronic structure theory.^ In particular, those states which arise from 
configurations which involve a doubly occupied d subshell have much larger 
correlation energies than those which do not.^ For example, the difference 
in correlation energies between the states arising from the d"s^ and 
occupancies is of the order of leV for the elements Sc-V, where the high-spin 
configuration never involves doubly occupied d orbitals, but increases 
to ~ 4eV for the elements Mn-Ni, in which the d""*"^ configuration has two 
additional filled d sublevels relative to <i"s^. The root of this problem lies in 
the very large Coulomb interaction and consequent correlation which occurs 
when two electrons are required to occupy the same, relatively small, 3d 
orbital. 

These errors can, in principle, be eliminated by configuration interaction 
(CI) calculations. Experience has shown,^"^''' however, that quite large one- 
electron basis sets including angular momenta at least through / functions are 
necessary in order to recover enough of the correlation energy to compute the 
energy differences reliably. Contracted basis sets of dimensions {10s7p4d3f) 
are necessary to approach the spdf limit. In addition, the many-electron 
basis must be quite large; CI limited to single and double substitutions yields 
errors of the order of 0.4-0. 7eV for the d^s^ — > d'^'^^s excitation energy for 
the elements to the right of the series (Mn-Cu), even when a large spdf basis 
is used and the results are corrected for relativistic contributions. 

The remaining error is associated with higher order excitations, and there- 
fore several size-extensive methods have been investigated. The most stan- 
dard approach, M0ller-Plesset perturbation theory, is not always adequate. 
In fact, for those atoms on the right hand side of the first transition series, the 
perturbation expansion fails dramatically and has not converged at fourth 
Q].(^g].8,i6,i7 (]VIP4). Coupled cluster techniques, in which certain classes of ex- 
citations are summed to infinite order, give the most reliable results to date. 
The Quadratic Configuration Interaction (QCI) approximation investigated 
by Raghavachari and Trucks^^'^^ reproduces the d"'s'^ d"+^s excitation en- 
ergies with a mean absolute error of 0.12cV at the QCISD(T) level of theory. 
A similar level of accuracy is obtained for the low-lying ionization potentials. 

A number of efforts have also been aimed at examining the performance of 
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density-functional theory (DFT).^^"^^ It is somewhat difficult to make a gen- 
eral statement concerning the level of accuracy achieved by these methods, 
because the magnitude of the errors in DFT approaches depend, as expected, 
on the specific exchange and correlation functionals used. A good review of 
this area, and the general status of DFT vs. Hartree-Fock schemes prior 
to 1987 is provided by Salahub.^ He points out that the DFT results also 
depend on "whether, and at what stage of the calculations spherical aver- 
aging is invoked." While most of the DFT results appear to be competitive 
with, or an improvement upon, Hartree-Fock theory, there are still rather 
large errors. Salahub concludes "the LSD [local spin density] method gen- 
erally presents a reasonable semi-quantitative picture and interprets trends 
correctly; however, it yields quantitative errors in relative energies as large 
as 1 eV or so." 

Nevertheless, a great deal of useful information has been obtained for 
inorganic complexes using the the local density approximation (LDA), i.e. 
the Slater exchange functional together with an electron correlation func- 
tional based on the properties of the homogeneous electron gas (e.g. the 
Vosko-Wilks-Nusair fit). Ziegler's recent review^^ is recommended, as are 
the articles in the book edited by Labanowski and Andzelm.^^ The most 
significant problem with this type of approach is the tendency of the LDA to 
overestimate molecular binding energies, sometimes by as much as 100%. A 
breakthrough in this regard has been the development of reliable "gradient- 
corrected" density functionals^^"^^ (sometimes referred to as non-local func- 
tionals). In particular, the gradient-corrected exchange functional of Becke^^ 
leads to much improved bond energies. 

Recently a number of studies of gradient-corrected DFT have begun to 
appear. ^^""^^ Geometries, vibrational frequencies, dipole moments, and other 
properties of small first row molecules are in reasonable agreement with ex- 
periment. A rule of thumb often mentioned is that these approaches are 
comparable in accuracy to MP2 theory. The gradient-corrected DFT ap- 
proaches have one clear advantage, however; the heats of atomization are 
remarkable both for their accuracy and for the fairly modest basis sets re- 
quired to achieve this accuracy. It seems apparent that the density converges 
much more rapidly with the one-electron basis than does the correlation en- 
ergy computed from Hartree-Fock based approaches. 

Given the success of gradient-corrected DFT for molecules composed of 
first row atoms, we were curious as to its performance for the excitation 



2 



energies and ionization potentials in the first transition series. The results 
of that study are reported in this paper. Section 2 reviews the methods 
used to study the Kohn-Sham equations, including issues regarding the basis 
sets, integration grids, and the details of an implementation of spin-restricted 
Kohn-Sham theory. The results are presented and discussed in Section 3, and 
the conclusions of this work reiterated in Section 4. 



The calculations described here were performed with the MESA suite of 
programs,^^ using a self-consistent Kohn-Sham (KS) procedure^^ with a fi- 
nite orbital (Cartesian-Gaussian) basis expansion For the closed shell species 
we have implemented the Kohn-Sham equations as described by Pople, Gill 
and Johnson.^'' Of particular note is the fact that this formulation leads to 
Fock-fike matrices which can be evaluated for gradient corrected functionals 
without evaluating the Hessian of the density. For the high-spin open-shell 
states we use a spin-restricted open-shell Kohn Sham (ROKS) procedure. 
Some form of this approach has apparently been used previously by Murray, 
Handy and Amos'^'^ to study the ^Bi state of CH2, but since they provide 
no details of the approach we outline our implementation briefiy below. We 
adopt the notation of Johnson et al.^^ 



One can view the Kohn-Sham equations as strictly analogous to the Hartree- 
Fock equations except for the replacement of the Hartree-Fock exchange op- 
erator with a local exchange-correlation potential. Given spin up and spin- 
down densities pa and pp, evaluated on a grid {vj} in space. 



where P"j, and P^^, are the familiar spin-up or spin-down density matrices and 
the (/)'s are the elements of the basis set, one can express a general functional 
of the density and its gradients as 



2 Methodology 



2.1 ROKS formalism 



(1) 



/ = /(Pq, P/3, 7aa, lap, 7/3/?) 



(2) 
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where the gradient invariants 7 are defined as 

7aa = I Vpal^ , laf) = ^ Pa ' Vp/J, 7/3/3 = | Vp/j]^ . (3) 

The exchange-correlation energy can be written as 

ExC = j fiPa,Pf3,laanal3, 1f30 )dT. (4) 

In order to solve the Kohn-Sham equations, one forms a Fock-like matrix as 



where 



^df ^ df ^ 



V{4>M dr, (6) 



and (i, j) G {(a, /?), (/?, a)}. The matrices /i and J are the usual one-electron 
and Coulomb matrices, respectively. In spin-unrestricted Kohn-Sham (UKS) 
calculations, the a and /? Fock matrices defined above are individually di- 
agonalized and the solutions iterated until self-consistency is achieved. For 
the spin-restricted high-spin open-shell calculations, we combine these oper- 
ators to obtain the analogue of a one-hamiltonian approach to Hartree-Fock 
theory. That is, if the orbitals are partitioned into closed, open, and virtual 
blocks, the matrix 

^0 Fco Fq \ 

Fo Fov (7) 
V Fo J 

is formed where 

Fo = ^{F^ + F^), 

Fco = Ff, (8) 
Fov = F". 

This operator is identical to Hamilton and Pulay's one-hamiltonian formula- 
tion of Hartree-Fock theory, except that the closed shell exchange matrix (K) 
is replaced by F^^^ and the open shell exchange matrix with F-^^"" — F-^'~^^ 
(the unpaired electrons are assigned alpha spin). The Kohn-Sham matrix 
equations are solved by the usual self-consistent techniques. The open-shell 
DIIS scheme of Hamilton and Pulay^^ was used to accelerate convergence. 
The Coulomb and one-electron terms are computed analytically. 
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2.2 Functionals 



For convenience our program calculates exchange and correlation functionals 
separately. For most of the present work we have used the gradient-corrected 
exchange functional of Becke^^ and the correlation functional of Lee, Yang 
and Parr.^^ We refer to this combination as B-LYP. For the copper atom we 
have also used Becke's functional with no correlation (B-nuU). 



2.2.1 Becke Exchange 

Becke's exchange functional, which was designed to improve upon the sim- 
ple local density approximation (LDA) by yielding the correct asymptotic 
behavior of the exchange energy, is given by 

fsiPa, Pl3, laa, 7/3/?) = pTqM + Pp^gixp) (9) 

with 

g^x) = ^f^Y' bx^ (10) 

2 \4ttJ 1 -|- 66a;sinh~^ x 

■^"^ ~ 4/3" \^^) 
Pa 

^/^ = (12) 

Note that the expressions for x^ and X/j given in Reference 37 contain a 
typographical error. The parameter h is given by Becke,^^ h = 0.0042. 



2.2.2 Lee- Yang-Parr correlation 

The Lee- Yang-Parr (LYP) correlation functional is derived from the correla- 
tion energy formula of CoUe and Salvetti,^^ derived from a consideration of 
short range effects in the two-particle density matrix. The functional itself, 
as transformed by Miehlich et al.,^^ is given by 

fLYp{Pa, Pp, laai 7«/3, 7/3/3) = 
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with 



LYP 



9^aa 
dfLYP 

dfLYP 

Hp) 



-2''/'^{3nr'abuip)p^p,ipT + pT) 

, dfiYP dfiYP , dfLYP 
+ — l(xa + -57 lap + ^ 7/3/9 



— —abuj{p) 



a/3 



^7/3/3 



Pa 



7TPaP/3 l-35(p)-[(5(p)-ll]^ 



-ahuj{p) {^PaP/3[47 - 7^(p)] - ^p^l 



-abui{p) 



P 



P/9 



l-3J(p)-[5(p)-ll]^ 



Cp-V3 + 



-11/3 



dp' 



-1/3 



1 + cZp-V3 • 



Pa 



(13) 



(14) 
(15) 
(16) 
(17) 
(18) 



The constants a,6, c, and d used were those used by Miehhch et al., ^ a — 
0.04918, h = 0.132, c = 0.2533, and d = 0.349. 



2.3 Grids and Integration Scheme 

The evaluation of the exchange- correlation contribution to the Fock matrix 
(Eq. 6) and the total energy (Eq. 4) is done by a numerical quadrature. 
In this paper we report atomic results but we have implemented the more 
general partitioning scheme of Becke.^^ For each atom in the molecule an 
atomic grid is generated consisting of concentric spherical shells centered 
on the atom. The shells are given radii according to the Euler-MacLaurin 
scheme described by Murray, et al.^^; each spherical shell has points and 
weights distributed on it according to the formula of Lebedev.^''' When the 
atomic grids are assembled into a molecular grid the weights are adjusted 
to partition the space into "fuzzy Voronoi polyhedra". For these atomic 
calculations the Lebedev angular grid of order 9 was used and is more than 
sufficient to integrate exactly all spherical harmonics which may appear in the 
solutions. The radial grid consisted of 100 points, and the "radial maximum" 
parameter chosen according to Slater's rules. In particular, we used values 
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of 4.59, 4.38, 4.20, 4.01, 3.82, 3.69, 3.53, 3.40, 3.27 and 3.16ao for the atoms 
Sc-Zn, respectively. The accuracy of the radial grid was tested on copper 
by comparing the results using both larger (150 and 200 radial points) and 
smaller (50 and 75 points) grids. The use of 100 radial points yields stability 
in the total energy (in a.u.) to four digits past the decimal point; increasing 
the grid to 150 points changes only the fifth digit, and using only 50 points 
changes the second digit. 

2.4 Basis Sets 

The Kohn-Sham orbitals were expanded in a primitive Cartesian Gaussian 
basis set of dimension (14s9p6d). This consists of Wachters' primitive set^^ 
augmented with the diffuse d function recommended by Hay.^^ All six 
Cartesian components of the d functions were retained. Since Wachters' 
primitive basis was optimized for Hartree-Fock calculations, and the detailed 
shape of the Kohn-Sham orbitals can, in principle, be different, some atten- 
tion was paid to the nature of the Kohn-Sham coefficients as an indication 
that the primitive basis should be modified. In general it appears that the 
Wachters/Hay set is adequate to expand the Kohn-Sham orbitals as well. In 
particular, most orbitals were described by two or more primitives with co- 
efficients of ~ 0.5. There were a few exceptions which might warrant further 
study. For example, whereas the Ti 2s Hartree-Fock orbital involves primi- 
tives 9 and 10 with coefficients of 0.5 and 0.6, the Ti 2s B-LYP orbital wants 
to be somewhat more diffuse as reflected by coefficients for primitives 9 and 
10 of 0.2 and 1.0. If one makes an even-tempered plot of the Ti s-space, 
it is apparent that there is a "gap", or missing s-function in the sequence 
which reflects the 2s /3s shell structure. An additional function would pre- 
sumably improve the description of this core orbital and reduce the total 
energy somewhat, but should have little effect on the valence electron distri- 
bution or relative energy differences. Additional examples of this sort occur 
for atoms on the right side of the series. 

We have also performed some initial investigations on contracted basis 
sets appropriate for use in the flrst transition series. These results are dis- 
cussed in Section 3.3. 
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2.5 Convergence issues 

For many of the cases described below we found that the simple ROKS 
procedure failed to converge. In most of these cases the cause of the failure 
was a tendency of the unpaired electrons to occupy different components of 
the 3d orbital from one iteration to the next. In most cases all that was 
required to avoid this was to employ symmetry constraints; for example, 
for the 3(P4s state of scandium we used the octahedral point group with 
inversion, and required that the two singly occupied d orbitals be those of the 
Eg irreducible representation, namely d2z2-x^-y2 and d^^-y^. An initial guess 
with the proper symmetry was formed by permuting the default initial guess 
vectors. At each step the Fock matrix was diagonalized in the symmetrized 
orbital basis, and the occupied orbitals chosen so that the number of occupied 
orbitals in each irreducible representation remained the same as for the initial 
guess. Enforcing Oh symmetry also prevents the s orbitals from mixing with 
the dz2 orbital. 

While using symmetry constraints was sufficient to ensure convergence 

for most states, there were a few stubborn cases, such as Sc'''(3(i4s). In these 
cases there were no symmetry groups which could both prevent mixing of 
the s and d orbitals and still prevent unpaired d electrons from moving be- 
tween degenerate symmetry-equivalent d orbitals. We resolved this problem 
by employing, for those cases where it was necessary, a maximum overlap 
condition in addition to the symmetry constraints: orbitals were filled ac- 
cording to the maximum overlap with the initial guess. The cases for which 
this was used were: Sc+(3fi4s), V(3dHs), Cr{3d^), Cr{3dHs^), Cr+(3(iMs), 
Fe{3dHs'^), Fe+{3dHs), and Co(3(i'^). 

3 Results 

3.1 Ionization and excitation energies 

The results of the B-LYP calculations of the ionization and excitation energies 
of the first transition series are given in Tables 1 and 2. The experimental 
values corrected for relativistic contributions are given in the final column 
of Tables 1 and 2. These values are simply the experimental results with an 
estimate of the differential relativistic energy subtracted out, as was done by 
Raghavachari and Trucks. The estimate is based on directly computed 
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relativistic contributions to the Hartree-Fock energies,^ scaled by an estimate 
for the effects of electron correlation. Figures 1-4 compare the B-LYP results 
with these corrected experimental values. 

Perusal of Tables 1 and 2 shows that the B-LYP approximation gives a 
generally reliable picture of the relative stability of the states of the first 
transition series. Consider first the d^s^ d"'^^s excitation energy shown 
in Figure 1. Note that the general "sawtooth" behavior in the excitation 
energy is faithfully reproduced by the restricted open-shell B-LYP method. 
Nevertheless, the ground state of the atom is predicted incorrectly in a few 
cases: vanadium, where B-LYP predicts the d'^s state to be more stable than 
d^s^; iron, where B-LYP yields d'^s lower than d^s"^; and cobalt, where B-LYP 
places d^s lower than d^s^. A general tendency of B-LYP to favor d"'~^^s over 
d^s configurations is apparent in Figure 1. The d"+^s states are generally 
predicted to be too stable by ~ 0.5eV. In V, Co and Fe the experimental 
splitting between (i"+^s and (i"s^ states is of this order or smaller, and so the 
bias leads to an incorrect ordering. This tendency to favor d^+^s states is 
common to the LDA as well. Harris and Jones^^ found a bias towards (i"+^s of 
~ leV. Theory and experiment are compared for the <i"s^ — > c?""*"^ excitation 
in Figure 2. Again, the qualitative features of the trends are reproduced 
well, but the configurations rich in li-electrons (or poor in s-electrons) are 
consistently predicted to be too stable, in this case by ~ 0.8eV. 

Figure 3 shows the results for the ionization potential d^s"^ —>■ d^s, a 
transition in which the number of d-electrons remains constant. Here, the 
B-LYP approximation is in much better agreement with experiment. This 
would suggest that the errors in the excitation energies of the neutral arise 
primarily from a tendency to overbind the d-electrons, and that the ioniza- 
tion potential of the s-electron is roughly correct. This conclusion is generally 
consistent with the results for the d^s^ — > ionization potential plotted 
in Figure 4. Most are underestimated consistent with a bias toward config- 
urations rich in d-electrons, although this rule of thumb would not predict 
the overestimates which occur at the far right of the series. 

The mean absolute errors in the B-LYP approximation for the various ex- 
citation energies are given in the final column of Table 3 , and the individual 
errors plotted in Figures 5 and 6. For purposes of comparison. Table 3 also 
contains the mean absolute errors as computed by Raghavachari and Trucks 
at the Hartree-Fock, M0ller-Plesset, and QCISD(T) levels of approximation. 
For the d^s^ d^s ionization potential, the mean absolute B-LYP error 
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is only 0.10 cV, and compares favorably with both the MP4(0.10eV) and 
QCISD(T)(0.09eV) results. The B-LYP error is somewhat larger (0.21eV) 
for the (i"s^ — > ionization, but is again competitive with the QCISD(T) 
approximation (0.14eV), and clearly superior to MP4. B-LYP does not per- 
form as well for the excitation energies of the neutral, with mean absolute 
errors of O.SleV and 0.76 for the d^Ks^ cf-^^s and d"'~^'^ transitions, 

respectively. The former can be compared with the QCISD(T) mean absolute 
error of 0.12eV. The overall mean B-LYP errors, it should be emphasized, 
are significantly smaller than the error from Hartrce-Fock calculations, which 
can be as large as 4eV for the (i"s^ (i""*"^ excitation. QCISD(T) results 
were not reported for this excitation. 

Thus far we have compared the gradient-corrected B-LYP DFT results 
with Hartree-Fock based approximations. It is also instructive to compare 
the B-LYP results with other gradient-corrected DFT functionals, such as 
the Becke-Perdew (BP) variant. Table 4 compares the mean absolute errors 
in B-LYP ionization potentials with the recent results of Zicglcr and Li,^° 
who examined both the Becke-Perdew approximation and the LDA (Slater 
exchange and Vosko-Wilk-Nusair correlation functionals). While it should 
be kept in mind that this is not a completely fair comparison of the func- 
tionals since Ziegler and Li used the spin-unrestricted version of Kohn-Sham 
theory and basis sets different from ours, it is interesting to note that both 
the B-LYP and BP functionals give quite similar results. It is also interesting 
that the LDA, while in worse agreement with experiment than either of the 
gradient-corrected approaches, performs rather well. In this series of ioniza- 
tion potentials, the gradient corrections do not appear to give dramatically 
improved results. 

This observation is also in accord with the recent work of Kutzler and 
Painter. For the s- ionization potentials, e.g., they find mean absolute er- 
rors of 0.22eV and 0.39eV for the gradient-corrected functionals of Langreth, 
Mchl and Hu (LMH)^^ and the generalized-gradient-approximation (GGA) 
of Perdew and Yue,^^ respectively. These can be compared with the B-LYP 
error of O.lOeV, the Becke-Perdew result of O.lGeV, and the LDA error of 
0.28eV. 

For the s — > d excitation energies in the neutral species, the LDA meann 
absolute error is ~ 0.85eV.^^ The B-LYP error found in the present work 
(O.SleV) is comparable to the GGA result (~ 0.6cV) of Kutzler and Painter, 
and significantly better than the LMH error (~ 1.2eV). 
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3.2 Analysis of B-LYP functional 

Although the B-LYP error is large for the d"'~^^s excitation energy, the 

fact that the qualitative features of the experimental trends are well repro- 
duced (cf. Figure 1) suggests that this might be traced to a systematic error 
which could be improved with future functionals. An immediate question, 
then, is where does the error reside, in the exchange or correlation functional? 
The origin of the error is not as easy to pin down as one might think at first. 
Consider the d^s^ d^^s transition in the copper atom. The difference in 
the exact exchange contribution between these two states can be extracted 
from Hartree-Fock calculations on the two. It turns out to be -5.62eV (see 
Table 5). The differential exchange energy from either the B-nuU or B-LYP 
calculation is very similar (-5.69 and -5.74eV, respectively). The difference 
in the correlation energy between the two states can be inferred from the 
Hartree-Fock calculations and the experimental results; it is -1.54eV. The 
correlation energy from the B-LYP calculations is only -0.16eV. One thus 
might expect a rather large total B-LYP error of -|-1.26eV (-0.1 2eV from the 
exchange and +1.38eV from the correlation energy error). The actual error 
in the B-LYP calculation for this case is -0.2eV. The reason for this is that 
the B-LYP one-electron and Coulomb contributions are different from the 
Hartree-Fock values. The self-consistency aspect of the calculation makes a 
direct examination of the error difficult. 

In order to make a more direct comparison, we also examined the results 
of the B-LYP functional being applied to the Hartree-Fock density. These 
numbers are also displayed in Table 5. The differential Becke exchange en- 
ergy (-6.99eV) is now quite different from the B-LYP value (-5.74eV) and 
in significant disagreement with the Hartree-Fock value (-5.62eV). The dif- 
ferential LYP correlation energy is hardly changed (-0.16 and -O.lSeV) and 
still significantly underestimated. Thus it appears that there are rather large 
errors in both functionals. The total B-LYP (HF) prediction is in rather good 
agreement with experiment, but this comes about because the overestimate 
of the exchange energy tends to cancel the underestimate of the correlation 
energy. Table 5 provides data for a similar analysis of the other states of 
interest in the copper atom. They are consistent with the conclusion reached 
above. Thus it appears that for Cu, the LYP functional underestimates 
the correlation energy, and that the Becke functional overestimates the ex- 
change energies. The total exchange-correlation energy is reproduced rather 
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accurately, and both the B-LYP and B-LYP(HF) approaches are in good 
agreement with experiment. 

The analysis and the discussion above leads one to conclude that the 
accuracy in the final B-LYP result arises from a rather fortuitous cancella- 
tion of errors in the exchange and correlation functionals. However, there 
appear to be some systcmatics to the error, as a study of Figures 5 and 6 
demonstrate. With the exception of the dP's^ — >■ dPs ionization potential, 
in which the error is roughly linear in n, the curves typically exhibit a saw- 
tooth behavior, in which the error at n = 5(Mn) suddenly becomes more 
severe. At first glance, one might attempt to associate the sudden change 
in the error for dP's^ — > dJ^^^s which occurs between Cr and Mn (77, = 4 and 
n = 5) to the sudden appearance of doubly occupied ci-subshells in the d^s 
state of Mn. For example, the error in the d'^s^ d^s excitation energy 
is ~ — 0.4cV in Cr, abruptly increasing to ~ — 0.9eV for the d^s'^ — > d^s 
excitation in Mn. Since the d^s state is predicted to be too stable relative 
to d^s^, one might conclude that the correlation functional overestimates the 
intra-pair d-d correlation energy. However, this argument would predict a 
similar increase in the error for the d^s'^ — > d^ excitation in Cr. Figure 5 
shows that this is not the case, and that the increase again occurs at Mn, 
this time in the d^s"^ d'^ excitation energy. Similarly, if one argues that 
the error is associated with the exchange functional and the "special" nature 
of the half-filled d^ shell, then one is hard-pressed to explain why a similar 
jump in the error is not apparent in the d^s^ — > d^ excitation energy of Ti. 
At present we do not understand this behavior and further work is clearly 
warranted for this problem. 

3.3 Basis Set Contractions 

In order to test the functionals. wc have attempted to eliminate many of the 
uncertainties associated with basis set incompleteness by using the fully un- 
contracted (14s9p6£i) basis of Wachters/Hay. For applications to molecules 
this basis is too large, and should be appropriately contracted. An obvious 
approach would be to use a general contraction scheme based on the Kohn- 
Sham orbitals obtained in this work. One might expect, however, that differ- 
ent contractions would be necessary for difi^erent variants of DFT; e.g., a set 
of contractions appropriate for LDA calculations, a different set for B-LYP, 
etc. For our initial investigation, we decided to test a contraction scheme 
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based on the Hartree-Fock orbitals. We therefore contracted the {14s9p6d) 
primitive basis into a [6s5p3d] basis using tlie general contraction scheme 
of Raffenetti^^ and the Hartree-Fock coefficients of Wachters.^^ Specifically, 
the inner parts of the Is, 2s, 2p, 3s, 3p and 3d orbitals were contracted via 
the HF coefficients, and the remaining, more diffuse primitives in each space 
were left free. These results are also presented in Tables 1 and 2. While 
the absolute energies were affected as expected, the excitation and ionization 
energies were generally changed by at most 0.02eV. The exceptions to this 
behavior occur for the Ni and Cu atoms, where some of the relative energies 
changed by as much as 0.12 eV. 

4 Conclusions 

We have examined the predictions of spin-restricted Kohn-Sham theory as 
regards the excitation and ionization energies of the members of the first 
transition series using the gradient-corrected density functionals of Becke 
and Lee, Yang and Parr. 

First of all, it is important to note that the qualitative features of the 
trends in excitation and ionization energies (Figs. 1-4) are faithfully repro- 
duced with a spin-restricted formulation of Kohn-Sham theory. It has some- 
times been suggested in the literature that a spin-unrestricted formulation is 
necessary to reproduce, e.g., the break in the d^s"^ — > d"'~^^s excitation energy 
which occurs at Cr (Fig. 1). This is clearly not the case. We believe this 
to be important, for a spin-restricted approach has the advantage that the 
solutions are eigenfunctions of spin and the uncontrolled spin contamination 
which freqently occurs in unrestricted Hartree-Fock (UHF) or Kohn-Sham 
(UKS) calculations on transition metal complexes is thereby avoided. We 
suspect that ROKS shares the disadvantages of ROHF: simple bonds will 
not always dissociate properly, molecules in which the qualitative electronic 
structure is best viewed as biradical in character may not be treated well, 
etc. In the present context, one might expect a significant spin polarization 
in the unrestricted formalism for those atomic states with a large number of 
unpaired d-electrons. While this can certainly lead to total energies which 
are significantly lower than the spin-restricted energies, it is interesting to 
note that the relative energies calculated by the ROKS procedure in this 
work are in close agreement with the UKS results of Ziegler and Li. A di- 
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rect comparison of ROKS vs. UKS is made difficult by the fact that Ziegler 
and Li used the Becke-Perdew functional whereas we examined Becke-LYP. 
The close agreement suggests, however, that neither the spin-polarization 
effects nor the different functionals used cause significant differences for the 
predictions of the ionization potentials of the first transition series. 

The quantitative agreement with experiment is also encouraging. Ioniza- 
tion energies are observed to be in good agreement with experiment, with 
mean absolute errors of O.lOeV for the dJ^s^ — > d^s ionization, and a mean 
absolute error of 0.21eV for (i"s^ — >■ d"'~^^. These errors are much smaller than 
those obtained by Hartrce-Fock, and compare favorably with the QCISD(T) 
results of Raghavachari and Trucks, which have mean absolute errors of 0.09 
and 0.14eV, respectively. The ROKS B-LYP errors in the excitation energies 
are larger, O.SleV for the — > d^'^^s excitation, as compared with an av- 
erage absolute error of 0.12eV from the QCISD(T) calculations. The B-LYP 
approximation, like other DFT variants, consistently places the d"'~^^s states 
too low compared with d^s"^. 

A comparison of the results using the primitive (14s9p6d) basis of Wachters 
versus a (6s5p3d) general contraction based on Hartree-Fock coefficients 
shows that the contracted basis is generally in excellent agreement with the 
fully uncontracted basis. This demonstrates that the description of the po- 
tential in the valence region of the atom due to the core electrons is essentially 
the same in B-LYP and Hartree-Fock calculations. This suggests that the 
relativistic effective core potentials developed for the Hartree-Fock problem 
may be applicable to the DFT methods as well. 

In summary, the spin-restricted Kohn-Sham calculations using the B-LYP 
functional look promising for calculations on transition metal complexes. The 
atomic errors are significantly smaller than the HF and MP approximations, 
particularly for those elements to the right of the row where the MP series is 
slow to converge. While not as accurate as the coupled-cluster techniques for 
this row, e.g. QCISD(T), B-LYP achieves a reasonable level of accuracy at 
a much reduced level of effort. Although it is clear that improved function- 
als are needed for quantitative accuracy, we expect the B-LYP functional, 
its variants or descendants to be increasingly used for electronic structure 
calculations on transition metal complexes. 
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Atom 


State 


Energy (au) 


(eV) 


Energy (au) 


(eV) 


Rel. Corr. 






(14s9p6c/) 


[6s5p3(i] 


Exp.16.17 


Sc 




-760.6369 


0.00 


-760.6310 


0.00 








-760.6132 


0.64 


-760.6071 


0.64 


1.33 






-760.5224 


3.12 


-760.5166 


3.11 


4.04 


Sc+ 


3d4s{^D) 
3d''{^F) 


-760.3992 


6.47 


-760.3932 


6.47 


6.54 




-760.3925 


6.65 


-760.3866 


6.65 


6.98 


Ti 


3dHs\^F) 


-849.3803 


0.00 


-849.3729 


0.00 






3dHs{^F) 
3d\^D) 


-849.3764 


0.11 


-849.3691 


0.10 


0.69 




-849.2851 


2.59 


-849.2781 


2.58 


3.17 


Ti+ 


3dHs{^F) 


-849.1326 


6.74 


-849.1252 


6.74 


6.80 




3d%^F) 


-849.1397 


6.55 


-849.1326 


6.54 


6.73 


V 


3dHs\^F) 


-943.9282 


0.00 


-943.9200 


0.00 






3dHs{^D) 


-943.9422 


-0.38 


-943.9341 


-0.38 


0.11 






-943.8724 


1.52 


-943.8646 


1.50 


2.24 


v+ 


3dHs{^F) 
3d\^D) 


-943.6711 


6.99 


-943.6628 


7.00 


7.03 




-943.6913 


6.45 


-943.6833 


6.44 


6.48 


Cr 


3dHs\^D) 


-1044.4196 


0.00 


-1044.4100 


0.00 






3dHsCS) 


-1044.4767 


-1.55 


-1044.4673 


-1.56 


-1.17 






-1011.3311 


2.10 


-1011.3213 


2.11 


3.14 


Cr+ 


3dHs{^D) 


-1044.1533 


7.25 


-1044.1442 


7.23 


7.24 






-1044.2144 


5.58 


-1044.2052 


5.57 


5.46 


Mn 




-1151.0244 


0.00 


-1151.0132 


0.00 






3dHs{^D) 


-1150.9863 


1.04 


-1150.9747 


1.05 


1.97 




3d''{*F) 


-1150.8728 


4.13 


-1150.8614 


4.13 


5.31 


Mn+ 


3dHs{'S) 


-1150.7514 


7.43 


-1150.7401 


7.43 


7.38 




3d%''D) 


-1150.7126 


8.48 


-1150.7007 


8.50 


8.92 


Fe 


3dHs\^D) 


-1263.7208 


0.00 


-1263.7145 


0.00 






3dUs{^F) 


-1263.7253 


-0.12 


-1263.7189 


-0.12 


0.65 




3d\^F) 


-1263.6219 


2.69 


-1263.6154 


2.70 


3.73 


Fe+ 


3dHs{'^D) 


-1263.4303 


7.90 


-1263.4240 


7.90 


7.84 




3d\^F) 


-1263.4434 


7.55 


-1263.4369 


7.55 


7.77 



Table 1: Results for B-LYP runs on the transition metals. 
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Atom State Energy (au) (eV) Energy (au) is.E (eV) Rel. Corr. 







{Us9p6d) 




[6s5p3d] 




Exp. 16-17 


Co 


2>dHs\^F) 


-1382.7912 


0.00 


-1382.7831 


0.00 






-1382.8076 


-0.45 


-1382.7996 


-0.45 


0.17 






-1382.7139 


2.10 


-1382.7057 


2.11 


2.95 


Co+ 


M\^F) 


-1382.4850 


8.33 


-1382.4769 


8.33 


8.20 




-1382.5179 


7.43 


-1382.5098 


7.43 


7.40 


Ni 


3dHs'('F) 


-1508.3439 


0.00 


-1508.3289 


0.00 






Sd^Asi^'D) 


-1508.3716 


-0.75 


-1508.3592 


-0.82 


-0.33 




Sd^XS) 


-1508.3199 


0.65 


-1508.3060 


0.62 


1.24 


Ni+ 


3dHs{^F) 


-1508.0225 


8.75 


-1508.0079 


8.73 


8.56 




3d\^D) 


-1508.0751 


7.31 


-1508.0612 


7.28 


7.08 


Cu 


3dHs\^D) 


-1640.5003 


0.00 


-1640.4797 


0.00 






Sd^Hs^S) 


-1640.5758 


-2.05 


-1640.5586 


-2.15 


-1.85 


Cu+ 


SdHsi^'D) 


-1640.1676 


9.05 


-1640.1444 


9.12 


8.92 




3d^%^S) 


-1640.2753 


6.12 


-1640.2590 


6.00 


5.65 


Zu 




-1779.1837 


0.00 


-1779.1656 


0.00 




Zn+ 


3d^Hs{^S) 
SdHs^^D) 


-1779.1345 


9.50 


-1779.1169 


9.49 


9.23 




-1778.8256 


17.91 


-1778.8048 


17.98 


17.85 



Table 2: Results for B-LYP runs on the transition metals (continued). 
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Ionization 


HF 


MP2 


MPS 


MP4 


QCISD(T) 


B-LYP 




1.43 


0.32 


0.31 


0.10 


0.09 


0.10 




0.75 


0.52 


0.43 


0.38 


0.14 


0.22 


Excitations 


dng2 _^ ^n+lg 


0.77 


0.49 


0.53 


0.41 


0.12 


0.51 


dng2 _^ 












0.76 



Table 3: Comparison of mean absolute deviations for the ionization and ex- 
citation energies of the first row transition metals. The B-LYP results are 
those of the current work, all others are taken from Raghavachari and Trucks 
(RT).^'''^^ Errors are relative to the experimental values with relativistic cor- 
rections. 



Ionization 


LDA'^ 


B-P" 


B-LYP'^ 




0.28 


0.16 


0.10 


d"s2 ^ ^n+l 


0.27 


0.23 


0.22 



Table 4: Comparison of mean absolute deviations (eV) for the ionization 
potentials of the first row transition metals. " Slater Exchange with VWN 
correlation functional^° ^ Becke Exchange with Perdew Correlation Func- 
tional. Present work. 



18 







HF 


B-null 


B-LYP 


B-LYP (HF) 






-5.62 


-5.69 


-5.74 


-6.99 


AE 

' ' ' 'corr 




-1.54 


0.00 


-0.16 


-0.18 


^Ecoul 




185.18 


153.30 


154.01 


185.18 






-179.88 


-149.50 


-150.16 


-179.88 


Etotal 






-1.89 


-2.05 


-1.86 




d^s 














3.40 


4.18 


4.23 


4.12 


AE 

' ' ' 'corr 




1.57 


0.00 


0.92 


0.92 


^Ecoul 




-255.57 


-259.55 


-264.33 


-255.57 


Ah 




259.52 


263.60 


268.23 


259.52 


Etotal 






8.23 


9.05 


8.98 














^Eg^ch 




-2.51 


-1.53 


-1.49 


2.90 


AE 

'—^^corr 




-0.45 


0.00 


0.40 


0.31 


^Ecoul 




-60.38 


-95.15 


-98.52 


-60.38 


Ahi 




69.00 


102.41 


105.73 


69.00 


Etotal 






5.73 


6.12 


6.02 



Table 5: Breakdown of energy contributions for B-LYP, B-null, HF and B- 
LYP(HF) (B-LYP using Hartree-Fock density) calculations on Cu. For the 
HF case Ecorr is the empirical correlation, Eg^p — Ehf 
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6 Figure Captions 



Energy (eV) 



-0.5 



-1.5 



-2.5 



B-LYP 
Experiment 




Figure 1: Plot of experimental and B-LYP values for the interconfigurational 

energy between the dP's^ and d^^^s states for the first transition series. Di- 
amonds are B-LYP values, plusses are experimental values with relativistic 
corrections. 
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5 

4.5 



3.5 



Energy (eV) 3 



1.5 



B-LYP 
Experiment 




Figure 2: Plot of experimental and B-LYP values for the interconfigurational 
energy between the dP's^ and d""^^ states for the first transition series. Di- 
amonds are B-LYP values, plusses are experimental values with relativistic 
corrections. 
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123456789 10 

n 

Figure 3: Plot of experimental and B-LYP values for the interconfigurational 
energy between the dT's^ and dP's states for the first transition series. Di- 
amonds are B-LYP values, plusses are experimental values with relativistic 
corrections. 
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Energy (eV) 



B-LYP 
Experiment 




Figure 4: Plot of experimental and B-LYP values for the interconfigurational 
energy between the dP's^ and dP'^^ states for the first transition series. Di- 
amonds are B-LYP values, plusses are experimental values with relativistic 
corrections. 
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Error (eV) 




123456789 



Figure 5: Plot of difference between experimental values and B-LYP values 
of excitation energies for the first row transition metals. Diamonds are the 
s^dP" — > sdP-^^ excitation energy, plusses are s^dP- — > 
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-0.5 




Figure 6: Plot of difference between experimental values and B-LYP values 
of ionization potentials for the first row transition metals. Diamonds are the 
s'^dP' — > sc?" ionization potential, plusses are s'^d^ — > d^'^^. 
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